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We investigate the role of migration patterns on the spread of epidemics in complex networks. We 
enhance the SIS-diffusion model on metapopulations to a nonlinear diffusion. Specihcally, individ¬ 
uals move randomly over the network but at a rate depending on the population of the departure 
patch. In the absence of epidemics, the migration-driven equilibrium is described by quantifying 
the total number of individuals living in heavily/lightly populated areas. Our analytical approach 
reveals that strengthening the migration from populous areas contains the infection at the early 
stage of the epidemic. Moreover, depending on the exponent of the nonlinear diffusion rate, epi¬ 
demic outbreaks do not always occur in the most populated areas as one might expect. 

PACS numbers: 89.75.Fb 87.23.Cc 


I. INTRODUCTION 

Density dependence is a crucial factor in dispersal 
movements of individuals in many species, be they an¬ 
imal or human. In turn, these movements are a key fac¬ 
tor in the fate of epidemic outbreaks occurring within 
spatially distributed populations [1]. 

Migration or dispersal, as the movement of individuals 
from one place to another, are affected by the presence 
of other conspecific members in two different ways. A 
positive density dependence reflects the fact that high 
population densities increase competition effects among 
their individuals and induce emigration from heavily pop¬ 
ulated locations. On the contrary, a negative density de¬ 
pendence corresponds to conspecific attraction, i.e., the 
tendency for individuals of a species to settle near con- 
specifics, and it has been claimed to be one of the reasons 
for leaving low-density areas [2] and for population ag¬ 
gregation patterns [3]. In general, there is empirical evi¬ 
dence that diffusion rates can both increase and decrease 
with population density. For instance, in some insect 
groups, crowded conditions lead to the appearance of a 
greater fraction of long-winged adults [4] and, in aphids, 
the appearance of winged individuals within populations 
of wingless adults [5]. On the other hand, negative den¬ 
sity dependence has been reported, for instance, in gulls, 
voles, and deers [2]. In humans, both positive and neg¬ 
ative dependences have been considered for prehistoric 
and historical human population dispersals [6]. 

Traditionally, populations are assumed to be con¬ 
tinuously distributed over space and, hence, nonlinear 
reaction-diffusion equations have been used to study the 
effect of the previous density dependencies on the spa¬ 
tial population dynamics [3, 7]. In this context, linear 
dependence of the diffusion coefficient on the population 
density has been proposed by several authors for animal 
dispersal [8, 9], whereas a power law with positive expo¬ 
nent has been used to describe the relationship between 
the insect dispersal rate and the population density [7]. 

Alternatively, the increasing fragmentation of habitats 


of many species, as well as the fact that, at a large scale, 
human travel can be described by flows among a set 
of discrete locations, make metapopulation models very 
useful for studying the dynamics of spatially subdivided 
populations and, also, for the analysis of spreading pro¬ 
cesses on top of these populations. Under this modeling 
approach, the spatial distribution of the whole popula¬ 
tion is described by a network of patches inhabited by 
local populations, as cities or habitats in a patchy land¬ 
scape, with migratory flows connecting them. 

The nature of these flows is, in fact, a key ingredient of 
metapopulation models. On the one hand, flows can be 
due to uniformly random migration, which assumes that 
any neighbouring patch of a “source” patch is reached 
with the same probability. In this case, the migration 
probability only depends on features of the source patch 
as, for instance, its population density or its degree or 
connectivity. On the other hand, migratory flows be¬ 
tween patches can depend on features of both source and 
destination patches. Such a non-uniformly random mi¬ 
gration is, in fact, assumed by the so-called gravity model 
of movements of people and goods traditionally used in 
social sciences [10, 11], and by the radiation model intro¬ 
duced in [12] to overcome some of its limitations. In an 
epidemic context, migration towards other patches can 
be even more related to the healthy conditions in the 
destination patches than to the conditions in the origin 
patch. In [13], for instance, a constant diffusion rate is 
assumed for each class of individuals, but the migration 
probability between two patches depends on the suscep¬ 
tible population in the destination patch. 

The goal of this paper is to deal with positive and 
negative density-dependent diffusion processes extending 
the equations for continuous-time epidemic dynamics on 
metapopulations derived in [14], which were based on the 
formalism of complex networks introduced in [15, 16]. In 
particular, we study the impact of density-dependent dif¬ 
fusion coefficients on the population distribution among 
heavily populated and lightly populated areas, how these 
distributions affect the epidemic growth, and the contri- 
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button of each local population to the epidemic spread¬ 
ing. 


II. THE MODEL 

The spatial arrangement of patches (areas) is described 
in terms of the connectivity (degree) distribution p{k) 
and the conditional probability P{k'\k) for a patch (node) 
of degree k to be connected to a patch of degree k'. More¬ 
over, within a patch of degree k, any individual has the 
same probability of leaving it throngh any of its k links, 
namely, 1/fc, which implies that the strength of the con¬ 
nections is independent of the travel distance between 
patches. 

According to these assumptions and denoting by 
PS,kit) and Pi, tit) the average number of suscepti¬ 
ble and infectious individuals in a patch of degree k 
(fc = 1, ..., /cmax), respectively, the following equations 
describe the epidemic dynamics over a metapopulation 

P's,kit) = (^p-I3cipk)^'j pi,k + Sipk-Ps,k) 

-Dsipk)ps,k+kY,DsiPk')Pik'\k)^ , 

k' 

< 

p'i,ki^) = (/5c(pfc)^ - pi,k - Spi,k 

-Diipk)pi,k+kY,DiiPk')Pik'\k)^ , 

,( 1 ) 

where /3 denotes the infection transmission probability 
through an infectious contact, p is the recovery rate, and 
equal birth and death rates 5 are assumed. The aver¬ 
age population size in a patch of degree k is Pkit) = 
PS,kit) + Pi, kit), and the contact rate c(pfc) is a non¬ 
decreasing density-dependent function, generalizing the 
two cases considered in [14, 17]: cip) = p (fully-mixed 
population) and c(p) = I (limited homogeneous mixing). 
Finally, the density-dependent diffusion rates of suscep¬ 
tible and infectious individuals are denoted by Dsipk) 
and Diipk), respectively. Note that it is natural to as¬ 
sume that the total outflow of individuals of each type 
in a patch, Dsip)p and Diip)p, are zero when p = 0 
(see Discussion for a comment on the implications of the 
violation of this hypothesis). 

The first term in each equation of (I) corresponds to 
the infection and recovery processes. The second term is 
the neutral demographic turnover. The last one corre¬ 
sponds to the migration/diffusion process which can be 
split into negative and positive terms. The former count 
the number of individuals leaving a patch of degree k, 
whereas the latter are the sum of the flows of individuals 
arriving at this patch from patches of degree k' , provided 
that patches of degree k are connected to patches of de¬ 
gree k' , i.e., for those k' such that P(fc'|fc) > 0. 

Once we know the solution of the system (1), the 
(expected) total number of susceptible and infected in¬ 
dividuals are Sit) = N ^j.pik) ps,kit) and lit) = 


^ '^kPi^) Pi,kit), respectively, where N is the number 
of nodes of the network. From Eq. (I) and assuming the 
consistency condition fcP(fc'|fc)p(fc) = fc'P(A:|fc')p(fc'), it 
follows that the total number of individuals is conserved 
in the metapopulation, i.e., S'(t) -I- /(t) = Np^, with p^ 
being the average number of individuals per patch. 

Denoting by ps and pi the vectors of the population 
distribution of each type in the metapopulation, system 
(1) can be written as 

p's = diag(/i - iScipk)^^) Pi 
Pk 

PiC-Id)<liagiDsipk))ps, 

< ( 2 ) 

p'l = diag(/3c(pfe)^^ - p) Pi 
Pk 

+ iC-Id)d\agiDiipk))pi, 

where p = p + S, C denotes the connectivity matrix 
with Ckk' = kPik'\k)/k', and diag(-) denotes a diago¬ 
nal matrix. C is non-negative and assumed to be irre¬ 
ducible. For the case of uncorrelated networks, Pik'\k) = 
k'pik')/{k) and, hence, Ckk' = kpik')/{k), where the 
brackets stand for the mean value. 


III. DIFFUSION WITHOUT EPIDEMICS 

In the absence of epidemic (p/ = 0), the system is 
driven by the diffusion process. To study the impact 
of diffusion rates on the population distribution pk, we 
consider the equation p's = iC — Id) diag(I?s(pfc)) ps and 
study the disease-free equilibrium pg- As C is irreducible, 
Vk = k, k = 1,..., fcmax 7 is the only positive eigenvector 
oiC—Id associated to the dominant eigenvalue A = 0. So, 
pI satisfies DsiPk)Pk = M k, with M being a suitable 
constant. Assuming that the total ontflow of individuals 
per patch P(p) := Dsip)p is continuous, strictly increas¬ 
ing, and satisfies the natural hypothesis P(0) = 0, the 
existence and uniqueness of a disease-free eqnilibrium 

pl=F-\Mk) (3) 

is guaranteed and, hence, M is computed from the nor¬ 
malizing condition ^kPi^)Pk ~ P'^■ worth noting 

that pj does not depend on the conditional probability 
P(fc'|fc), so it is independent of the degree-degree corre¬ 
lations, and, moreover, that it increases with k since the 
total outflow per patch A(p) is an increasing function. 

To deal with both positive and negative density depen¬ 
dencies, we will assume the following typical nonlinear 
form for the diffusion rate Dg (see [3] and the references 
therein) 

Dsip) = (^) , a > -1, (4) 

with Dg > 0 (units of time“^), which guarantees the 
previous hypotheses on the total outflow per patch 
F(p). Regarding to the dimensionless parameter a, for 
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that a patch of connectivity k is heavily populated (HP) 
if its population is greater than the average number of 
individuals per patch in the metapopulation, i.e., when 
Pk > Otherwise, a patch is said to be lightly popu¬ 
lated {pk < p°)- 

When the disease-free equilibrium is computed from 
(3), HP patches turn out to be those with connectivity 
larger than k* := Ds{p^)p^/M, with M the normalizing 
constant for the profile. To illustrate heavily/lightly pop¬ 
ulated scenarios, we assume the form (4) for the diffusion 
rate Ds and a scale-free network with degree distribution 
p{k) = k > fcinin, 7 > 2. From (5), HP 

patches are those with connectivity 

^ ^ A: ■ 

I '^min 

(jO — 1 J 


k>k* = 


FIG. 1. (Color online) Percentage of population in heavily 
populated (HP) patches, ((ta — l)/a;)““^ x 100%, with uj = 
(7 — 1)(1-|-q:) > 1, as a function of the dimensionless migration 
exponent a in (4), and degree distribution p{k) ~ k ' with 
fcmin = 3 and 7 = 3 . Solid dots correspond to the profiles for 
a = —1/3 (heavily populated scenario), q = 0 when the ratio 
of the individuals living in heavily/lightly populated patches 
is one, and a = 1 (lightly populated scenario). 


— 1 < a < 0 we have the scenario in which the propen¬ 
sity to emigrate is higher in lightly populated patches 
{p < p°). On the other hand, an exponent a > 0 mod¬ 
els the scenario in which the propensity to emigrate is 
higher in heavily populated patches (p > p°), and it has 
been used for modelling animal dispersal in [7] (Section 
11.3). Finally, a = 0 recovers the constant diffusion rate 
considered in [14, 17]. 

Taking (4), the expression (3) for the disease-free equi¬ 
librium reads 


Pk Q,l/(l+a)\iP ■ 


(5) 


Note that this power-law distribution arises indepen¬ 
dently of the network topology. Moreover, other types 
of profiles are possible by just changing the form in (4). 
It is worth mentioning that, in terms of our notation, 
the same profile as (5) is found in [18] by assuming a 
constant overall diffusion rate from a patch, but dif¬ 
fusion rates between pairs of patches depending on the 
degrees of the involved nodes. Therefore, more than one 
mechanism of dispersal can explain power-law profiles of 
pI (see Discussion). 

Beyond the richer variety of population profiles found 
so far, it is interesting to know to what extent these sta¬ 
tionary profiles correspond to metapopulations whose in¬ 
dividuals live mostly in heavily populated patches or, on 
the contrary, in lightly populated ones. Precisely, we say 


where w = ( 7 —l)(l-|-a), and the total population in HP 
patches, in the limit of very large networks, is 

nJ pik)pldk = Np° , (6) 

with w > 1 to assure the convergence of the integrals. 
In the situation above and for w = 2, the percentages of 
individuals living in both types of patches are equal to 
50%. That would correspond, e.g., to the case of constant 
diffusion rate a = 0 and 7 = 3 as in [17]. However, other 
migration patterns allow us to describe a wider spectrum 
of population distributions, specifically, those with a per¬ 
centage of population living in HP patches higher (w < 2) 
or lower (w > 2) than 50%, see Fig. I for the case 7 = 3 
and varying a. For any exponent 7 > 2 we get analo¬ 
gous pictures as in Fig. 1. In particular, from ( 6 ) and 
for a fixed exponent 7 > 2 of the degree distribution, the 
migration exponent a can be used as a tuning parame¬ 
ter to shape the profile to a metapopulation with a given 
percentage of population in HP patches. 

However, the percentage of population living in HP 
patches, namely, ((w — x 100%, always lies in 

(36.79%, 100%). This is due to the fact that this percent¬ 
age is decreasing in w, with the limit values 100 e“^% 
and 100% obtained when w —>■ 00 and w —)• 1 , respec¬ 
tively. Notice that it is possible to obtain 100% of the 
total population in highly populated patches because the 
computation in ( 6 ) assumes infinite network sizes. 


IV. EARLY STAGE OF THE EPIDEMIC 

Now, let us focus on the effect of the migration patterns 
on the onset of the epidemic. The initial epidemic growth 
rate is given by the dominant eigenvalue (i.e. spectral 
bound) Ai of the Jacobian matrix of (2) at the disease- 
free equilibrium. This matrix can be written in blocks 
as: 
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migration exponent a 



FIG. 2. (Color online) Comparison of the initial epidemic growth rate (time”*^) for two contact patterns: c(p) = 100p/p° 
(left) and c(p) = 1000p/(p° + p) (right). Growth rate computed as its actual value (dots) given by the dominant eigenvalue 
Ai of the Jacobian matrix (fcmin = 3, fcmax = 220), and its estimation (solid line) given by left-hand side of (8). Results 
obtained for uncorrelated scale-free networks {-y = 3, N = 5000 nodes) and parameters: /3 = 0.015 (dimensionless); p = 1, 
Ds{p) = Di{p) — 0.5(p/p°)“ (units of time“^), and p° — 100 individuals. 


/(C-/d)-diag(F'(p*)) 

0 (C - Id) 


where the total outflow per patch F{p) = Ds{p)p is as¬ 
sumed to be differentiable. An exponential initial growth 
of the infectious population occurs when the disease-free 
equilibrium becomes unstable, i.e., when Ai > 0. The 
block structure of the Jacobian matrix and the fact that 
the dominant eigenvalue of the first block is 0, imply that 
Ai is given by the dominant eigenvalue of the fourth block 
when the latter is positive. 

For a constant contact rate, c(p) = cq, the dominant 
eigenvalue is equal to Ai = max{0, /3co — p} and it turns 
out to be independent of the migration pattern. On the 
contrary, when the contact rate is density dependent, a 
sufficient condition for Ai > 0 to be fulfilled is that 

max{^c(pfc) - p - (1 - P{k\k))Di{pl)} > 0 (8) 

k 

which follows after grouping all diagonal terms in the 
fourth block (see appendix A.8 in [19] for details). So, 
left-hand side of (8) is a lower bound of Ai since it ne¬ 
glects off-diagonal terms of the connectivity matrix. 

To have a global description of the role of the migra¬ 
tion, we have computed the initial growth of the epidemic 
Ai as the largest eigenvalue of the fourth block of the Ja¬ 
cobian matrix, for different migration exponents a and 
we have checked the goodness of its estimation given by 
the left-hand side of (8). Fig. 2 shows the fit for two con¬ 
tact patterns, namely, nonlimited and saturated number 
of contacts per unit of time. The first is the standard 
mass action c(p) = cop and in the second we consider a 


-diag (^c(p^)-p) \ 

diag {Di{pl)) + diag (/3c(p^) - p) y ’ 


saturated contact rate of the form c(p) = cop/(p° + p), 
Co > 0, which comes from the assumption that the time 
an individual is available for contacts is limited [19]. In 
both cases shown in Fig. 2, the initial growth rate of the 
epidemic decreases as the migration exponent increases. 

Recalling the definition of the basic reproduction num¬ 
ber of a local population as the number of secondary 
infections produced by a typical individual in a wholly 
susceptible population, and computed as infection prob¬ 
ability X contact rate x average infectious period, we can 
rewrite (8) as 


max 

k 


_/3c(Pfc)_ 

p + [I - P{k\k))Di{pl) 


> 1 . 


(9) 


So, for a given k, the ratio in (9) is a lower bound of the 
basic reproduction number of the populations living in 
patches with connectivity k, since the immigration from 
patches with connectivities k' ^ k\s neglected. Note that 
the factor 1—P{k\k) = Pik'\k) in the denominator 

accounts for those individuals who emigrate to patches 
with other connectivities. Therefore, local epidemic out¬ 
breaks will undoubtedly take place in those patches with 
connectivity k where the ratio in (9) is larger than one. 

In contrast to the case of constant diffusion [17], for 
density-dependent contact rates c(p) and for Dj(p) = 
tDs{p), t > 0, with Dg of the form (4), numerical com¬ 
putations show that the maximum in (9) is not necessar¬ 
ily attained at populations with the largest connectivity. 
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FIG. 3. (Color online) Dimensionless ratio in (9) for uncorrelated scale-free networks (7 = 3 and femin = 3) and diffusion 
Ds{p) = Di{p) — ■ {p/p^)°‘- Left panel, from left to right: migration exponent a = I, 1.7, and 3. Right panel: filled contour 

plot of the ratio in (9) showing that the maximum value decreases as a increases for the considered range of the parameter 
values. Parameters: /3 = 0.015 (dimensionless); p ~ 1, c{p) = 100p/p°, D° = 0.5 (units of time“^), and p° = 100 individuals. 


In some cases, this maximum is attained at the minimum 
degree feminj and, in others, at an intermediate value of k 
(see left panel in Fig. 3 for an illustration of these cases). 
In addition, the right panel of Fig. 3 is a plot of the 
ratio in (9) as a function of both the degree k and the 
migration exponent a, and shows that the higher the mi¬ 
gration exponent a, the lower the maximum value (9), 
for the considered range of parameters values. 


V. DISCUSSION 

Using the modelling framework introduced in [14], a 
density-dependent diffusion rate D{p) has been consid¬ 
ered to model positive and negative density dependencies 
of individuals’ migratory movements within a metapop¬ 
ulation [2]. In order to analyse the effect of these de¬ 
pendencies on the dispersal process, D{p) has been as¬ 
sumed to be a power-law function of p with an exponent 
a > —1. This lower bound for the allowed values of 
a guarantees both that the flow of individuals leaving 
a patch, D{p)p, strictly increases with its population p, 
and that Zl(p)p 0 as p —^ 0 (nobody migrates from 
where nobody lives). 

Under these two natural hypotheses on D{p)p, the 
disease-free equilibrium is determined by the diffusion 
process (which does not need to be necessarily described 
in terms of a power-law diffusion rate) and it is computed 
from the eigenvector associated to the largest eigenvalue 
of the connectivity matrix. With this respect, it is in¬ 
teresting to observe that if these hypotheses are violated 
as by assuming that the number of individuals traveling 
out of a patch, D{p)p, is independent from its population 
size p, the diffusion process does not determine the pro¬ 


file of the disease-free equilibrium. Such an assumption, 
indeed, leads to a diffusion rate D inversely proportional 
to the patch population p, and to an equilibrium which 
is uniquely determined by the initial distribution Pfc(O) 
(see Sect. 3.2 in [18]). 

On the other hand, when D[p) ^ p“ (a > —1), the 
disease-free equilibrium is pj ^ fci/(i+“) independently 
of the network topology. So, a power-law profile for p^ 
has been obtained under the only assumption of nonlin¬ 
ear diffusion. Remarkably, the same type of profile has 
been obtained for uncorrelated networks in [18] under a 
different physical principle, namely, assuming that the 
migratory flow ojkk' between an origin patch of degree k 
and a destination patch of degree k' is Wkk' = wo{kk'y, 
6 > 0. From these flows, the diffusion rate between nodes 
of degree k and k' is taken to be equal to D°Wkk'/Tk with 
Tk = kJ2k' Pik'\k)wkk', which ensures a constant diffu¬ 
sion rate from any patch in the metapopulation [18]. 

Which of these physical assumptions is mainly at work 
in a random dispersal process depends on the species we 
are interested in. In general, it is expectable that many 
insects and other animal species experience only local 
environmental conditions and, hence, density dependen¬ 
cies at the origin location would be the main factor in 
random dispersal within a metapopulation. On the con¬ 
trary, humans migrate with knowledge of the political, 
economic and environmental conditions in the potential 
arrival patches. Under this circumstance, it is suitable 
to assume that the migratory flow between two patches 
depends on features of the involved patches, and, in fact, 
it is what is assumed in many human mobility models 
[11, 12]. So, distinct physical mechanisms can be behind 
the observed scaling laws in the topological features of 
mobility networks [10]. 
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A way to quantify the effect of different density depen¬ 
dencies is by computing the percentage of individuals liv¬ 
ing in heavily populated patches (i.e., above the average 
population per patch). This aggregate measure of the 
population over the whole metapopulation is a comple¬ 
mentary description to that given by In particular, 
for infinite scale-free networks and assuming D{p) ^ 
this percentage turns out to be bounded from below by 
36.79%, which follows as a limit case of the set of power- 
law distributions, e.g. when the migration exponent a 
is large enough, i.e., when there is a high propensity to 
emigrate from heavily populated patches. 

From an epidemiological point of view, the main nov¬ 
elty of our results is to show that migration patterns de¬ 
termine where epidemic outbreaks can take place with 
higher probability. In particular, for large enough mi¬ 
gration exponents, the onset of an epidemic may not be 
triggered by infectious individuals living in large size pop¬ 
ulations, as it could be expected from the fact that 


increases with the connectivity of the sites. This hap¬ 
pens when individuals, also the infectious ones, have a 
high propensity to move from heavily populated areas. 
Such a higher number of infectious individuals arriving 
at patches with lower population sizes makes more likely 
the occurrence of an outbreak in these patches. But, at 
the same time, their lower sizes lead to lower values of the 
basic reproduction number because the contact rate c(p) 
increases with p (see Fig. 3). Therefore, in contrast to 
what is usually claimed about the role of the long-range 
mixing in exacerbating epidemics [20], we have seen that 
strengthening the migration from heavily populated ar¬ 
eas can help to contain an epidemic by causing its ap¬ 
pearance in less populated areas, for which the value of 
the basic reproduction number is significantly smaller. 
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